home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / INTERNET / SITES / LITTLE / P3SRC.ZIP / ATARI / POLY.C < prev    next >
Encoding:
C/C++ Source or Header  |  1997-01-18  |  29.5 KB  |  1,575 lines
  1. /****************************************************************************
  2. *                poly.c
  3. *
  4. *  This module implements the code for general 3 variable polynomial shapes
  5. *
  6. *  This file was written by Alexander Enzmann.  He wrote the code for
  7. *  4th - 6th order shapes and generously provided us these enhancements.
  8. *
  9. *  from Persistence of Vision(tm) Ray Tracer
  10. *  Copyright 1996 Persistence of Vision Team
  11. *---------------------------------------------------------------------------
  12. *  NOTICE: This source code file is provided so that users may experiment
  13. *  with enhancements to POV-Ray and to port the software to platforms other 
  14. *  than those supported by the POV-Ray Team.  There are strict rules under
  15. *  which you are permitted to use this file.  The rules are in the file
  16. *  named POVLEGAL.DOC which should be distributed with this file. If 
  17. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  18. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  19. *  Forum.  The latest version of POV-Ray may be found there as well.
  20. *
  21. * This program is based on the popular DKB raytracer version 2.12.
  22. * DKBTrace was originally written by David K. Buck.
  23. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  24. *
  25. *****************************************************************************/
  26.  
  27. #include "frame.h"
  28. #include "vector.h"
  29. #include "povproto.h"
  30. #include "bbox.h"
  31. #include "polysolv.h"
  32. #include "matrices.h"
  33. #include "objects.h"
  34. #include "poly.h"
  35. #include "povray.h"
  36.  
  37. /*
  38.  * Basic form of a quartic equation:
  39.  *
  40.  *  a00*x^4    + a01*x^3*y   + a02*x^3*z   + a03*x^3     + a04*x^2*y^2 +
  41.  *  a05*x^2*y*z+ a06*x^2*y   + a07*x^2*z^2 + a08*x^2*z   + a09*x^2     +
  42.  *  a10*x*y^3  + a11*x*y^2*z + a12*x*y^2   + a13*x*y*z^2 + a14*x*y*z   +
  43.  *  a15*x*y    + a16*x*z^3   + a17*x*z^2   + a18*x*z     + a19*x       + a20*y^4   +
  44.  *  a21*y^3*z  + a22*y^3     + a23*y^2*z^2 + a24*y^2*z   + a25*y^2     + a26*y*z^3 +
  45.  *  a27*y*z^2  + a28*y*z     + a29*y       + a30*z^4     + a31*z^3     + a32*z^2   + a33*z + a34
  46.  *
  47.  */
  48.  
  49.  
  50.  
  51. /*****************************************************************************
  52. * Local preprocessor defines
  53. ******************************************************************************/
  54.  
  55. #define DEPTH_TOLERANCE 1.0e-4
  56. #define INSIDE_TOLERANCE 1.0e-4
  57. #define ROOT_TOLERANCE 1.0e-4
  58. #define COEFF_LIMIT 1.0e-20
  59. #define BINOMSIZE 40
  60.  
  61.  
  62.  
  63. /*****************************************************************************
  64. * Local typedefs
  65. ******************************************************************************/
  66.  
  67.  
  68.  
  69. /*****************************************************************************
  70. * Static functions
  71. ******************************************************************************/
  72.  
  73. static int intersect PARAMS((RAY *Ray, int Order, DBL *Coeffs, int Sturm_Flag,
  74. DBL *Depths));
  75. static void normal0 PARAMS((VECTOR Result, int Order, DBL *Coeffs,
  76. VECTOR IPoint));
  77. static void normal1 PARAMS((VECTOR Result, int Order, DBL *Coeffs,
  78. VECTOR IPoint));
  79. static DBL inside PARAMS((VECTOR IPoint, int Order, DBL *Coeffs));
  80. static int intersect_linear PARAMS((RAY *ray, DBL *Coeffs, DBL *Depths));
  81. static int intersect_quadratic PARAMS((RAY *ray, DBL *Coeffs, DBL *Depths));
  82. static int factor_out PARAMS((int n, int i, int *c, int *s));
  83. static long binomial PARAMS((int n, int r));
  84. static void factor1 PARAMS((int n, int *c, int *s));
  85.  
  86. /* unused
  87. static DBL evaluate_linear PARAMS((VECTOR P, DBL *a));
  88. static DBL evaluate_quadratic PARAMS((VECTOR P, DBL *a));
  89. */
  90.  
  91. static int All_Poly_Intersections PARAMS((OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack));
  92. static int Inside_Poly PARAMS((VECTOR IPoint, OBJECT *Object));
  93. static void Poly_Normal PARAMS((VECTOR Result, OBJECT *Object, INTERSECTION *Inter));
  94. static void *Copy_Poly PARAMS((OBJECT *Object));
  95. static void Translate_Poly PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  96. static void Rotate_Poly PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  97. static void Scale_Poly PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  98. static void Transform_Poly PARAMS((OBJECT *Object, TRANSFORM *Trans));
  99. static void Invert_Poly PARAMS((OBJECT *Object));
  100. static void Destroy_Poly PARAMS((OBJECT *Object));
  101.  
  102. /*****************************************************************************
  103. * Local variables
  104. ******************************************************************************/
  105.  
  106. METHODS Poly_Methods =
  107. {
  108.   All_Poly_Intersections,
  109.   Inside_Poly, Poly_Normal, Copy_Poly,
  110.   Translate_Poly, Rotate_Poly,
  111.   Scale_Poly, Transform_Poly, Invert_Poly, Destroy_Poly
  112. };
  113.  
  114.  
  115.  
  116. /* The following table contains the binomial coefficients up to 15 */
  117.  
  118. static int binomials[15][15] =
  119. {
  120.   {1,  0,  0,  0,   0,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0},
  121.   {1,  1,  0,  0,   0,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0},
  122.   {1,  2,  1,  0,   0,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0},
  123.   {1,  3,  3,  1,   0,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0},
  124.   {1,  4,  6,  4,   1,   0,   0,   0,   0,   0,   0,  0,  0,  0,  0},
  125.   {1,  5, 10, 10,   5,   1,   0,   0,   0,   0,   0,  0,  0,  0,  0},
  126.   {1,  6, 15, 20,  15,   6,   1,   0,   0,   0,   0,  0,  0,  0,  0},
  127.   {1,  7, 21, 35,  35,  21,   7,   1,   0,   0,   0,  0,  0,  0,  0},
  128.   {1,  8, 28, 56,  70,  56,  28,   8,   1,   0,   0,  0,  0,  0,  0},
  129.   {1,  9, 36, 84, 126, 126,  84,  36,   9,   1,   0,  0,  0,  0,  0},
  130.   {1, 10, 45,120, 210, 252, 210, 120,  45,  10,   1,  0,  0,  0,  0},
  131.   {1, 11, 55,165, 330, 462, 462, 330, 165,  55,  11,  1,  0,  0,  0},
  132.   {1, 12, 66,220, 495, 792, 924, 792, 495, 220,  66, 12,  1,  0,  0},
  133.   {1, 13, 78,286, 715,1287,1716,1716,1287, 715, 286, 78, 13,  1,  0},
  134.   {1, 14, 91,364,1001,2002,3003,3432,3003,2002,1001,364, 91, 14,  1}
  135. };
  136.  
  137. static DBL eqn_v[3][MAX_ORDER+1], eqn_vt[3][MAX_ORDER+1];
  138.  
  139.  
  140.  
  141.  
  142. /*****************************************************************************
  143. *
  144. * FUNCTION
  145. *
  146. *   All_Poly_Intersections
  147. *
  148. * INPUT
  149. *   
  150. * OUTPUT
  151. *   
  152. * RETURNS
  153. *   
  154. * AUTHOR
  155. *
  156. *   Alexander Enzmann
  157. *   
  158. * DESCRIPTION
  159. *
  160. *   -
  161. *
  162. * CHANGES
  163. *
  164. *   -
  165. *
  166. ******************************************************************************/
  167.  
  168. static int All_Poly_Intersections(Object, Ray, Depth_Stack)
  169. OBJECT *Object;
  170. RAY *Ray;
  171. ISTACK *Depth_Stack;
  172. {
  173.   POLY *Poly = (POLY *) Object;
  174.   DBL Depths[MAX_ORDER], len;
  175.   VECTOR  IPoint;
  176.   int cnt, i, j, Intersection_Found, same_root;
  177.   RAY New_Ray;
  178.  
  179.   /* Transform the ray into the polynomial's space */
  180.  
  181.   MInvTransPoint(New_Ray.Initial, Ray->Initial, Poly->Trans);
  182.   MInvTransDirection(New_Ray.Direction, Ray->Direction, Poly->Trans);
  183.  
  184.   VLength(len, New_Ray.Direction);
  185.   VInverseScaleEq(New_Ray.Direction, len);
  186.  
  187.   Intersection_Found = FALSE;
  188.  
  189.   Increase_Counter(stats[Ray_Poly_Tests]);
  190.  
  191.   switch (Poly->Order)
  192.   {
  193.     case 1:
  194.  
  195.       cnt = intersect_linear(&New_Ray, Poly->Coeffs, Depths);
  196.  
  197.       break;
  198.  
  199.     case 2:
  200.  
  201.       cnt = intersect_quadratic(&New_Ray, Poly->Coeffs, Depths);
  202.  
  203.       break;
  204.  
  205.     default:
  206.  
  207.       cnt = intersect(&New_Ray, Poly->Order, Poly->Coeffs, Test_Flag(Poly, STURM_FLAG), Depths);
  208.   }
  209.  
  210.   if (cnt > 0)
  211.   {
  212.     Increase_Counter(stats[Ray_Poly_Tests_Succeeded]);
  213.   }
  214.  
  215.   for (i = 0; i < cnt; i++)
  216.   {
  217.     if (Depths[i] > DEPTH_TOLERANCE)
  218.     {
  219.       same_root = FALSE;
  220.  
  221.       for (j = 0; j < i; j++)
  222.       {
  223.         if (Depths[i] == Depths[j])
  224.         {
  225.           same_root = TRUE;
  226.  
  227.           break;
  228.         }
  229.       }
  230.  
  231.       if (!same_root)
  232.       {
  233.         VEvaluateRay(IPoint, New_Ray.Initial, Depths[i], New_Ray.Direction);
  234.  
  235.         /* Transform the point into world space */
  236.  
  237.         MTransPoint(IPoint, IPoint, Poly->Trans);
  238.  
  239.         if (Point_In_Clip(IPoint, Object->Clip))
  240.         {
  241.           push_entry(Depths[i] / len,IPoint,Object,Depth_Stack);
  242.  
  243.           Intersection_Found = TRUE;
  244.         }
  245.       }
  246.     }
  247.   }
  248.  
  249.   return (Intersection_Found);
  250. }
  251.  
  252.  
  253.  
  254. /*****************************************************************************
  255. *
  256. * FUNCTION
  257. *
  258. *   evaluate_linear
  259. *
  260. * INPUT
  261. *   
  262. * OUTPUT
  263. *   
  264. * RETURNS
  265. *   
  266. * AUTHOR
  267. *
  268. *   Alexander Enzmann
  269. *   
  270. * DESCRIPTION
  271. *
  272. *   -
  273. *
  274. * CHANGES
  275. *
  276. *   -
  277. *
  278. ******************************************************************************/
  279.  
  280. /* For speedup of low order polynomials, expand out the terms
  281.    involved in evaluating the poly. */
  282. /* unused
  283. static DBL evaluate_linear(P, a)
  284. VECTOR P;
  285. DBL *a;
  286. {
  287.   return(a[0] * P[X]) + (a[1] * P[Y]) + (a[2] * P[Z]) + a[3];
  288. }
  289. */
  290.  
  291.  
  292.  
  293. /*****************************************************************************
  294. *
  295. * FUNCTION
  296. *
  297. *   evaluate_quadratic
  298. *
  299. * INPUT
  300. *   
  301. * OUTPUT
  302. *   
  303. * RETURNS
  304. *   
  305. * AUTHOR
  306. *
  307. *   Alexander Enzmann
  308. *   
  309. * DESCRIPTION
  310. *
  311. *   -
  312. *
  313. * CHANGES
  314. *
  315. *   -
  316. *
  317. ******************************************************************************/
  318.  
  319. /*
  320. static DBL evaluate_quadratic(P, a)
  321. VECTOR P;
  322. DBL *a;
  323. {
  324.   DBL x, y, z;
  325.  
  326.   x = P[X];
  327.   y = P[Y];
  328.   z = P[Z];
  329.  
  330.   return(a[0] * x * x + a[1] * x * y + a[2] * x * z +
  331.          a[3] * x     + a[4] * y * y + a[5] * y * z +
  332.          a[6] * y     + a[7] * z * z + a[8] * z     + a[9]);
  333. }
  334. */
  335.  
  336.  
  337.  
  338. /*****************************************************************************
  339. *
  340. * FUNCTION
  341. *
  342. *   factor_out
  343. *
  344. * INPUT
  345. *   
  346. * OUTPUT
  347. *   
  348. * RETURNS
  349. *   
  350. * AUTHOR
  351. *
  352. *   Alexander Enzmann
  353. *   
  354. * DESCRIPTION
  355. *
  356. *   Remove all factors of i from n.
  357. *
  358. * CHANGES
  359. *
  360. *   -
  361. *
  362. ******************************************************************************/
  363.  
  364. static int factor_out(n, i, c, s)
  365. int n, i, *c, *s;
  366. {
  367.   while (!(n % i))
  368.   {
  369.     n /= i;
  370.  
  371.     s[(*c)++] = i;
  372.   }
  373.  
  374.   return(n);
  375. }
  376.  
  377.  
  378.  
  379. /*****************************************************************************
  380. *
  381. * FUNCTION
  382. *
  383. *   factor1
  384. *
  385. * INPUT
  386. *   
  387. * OUTPUT
  388. *   
  389. * RETURNS
  390. *   
  391. * AUTHOR
  392. *
  393. *   Alexander Enzmann
  394. *   
  395. * DESCRIPTION
  396. *
  397. *   Find all prime factors of n. (Note that n must be less than 2^15.
  398. *
  399. * CHANGES
  400. *
  401. *   -
  402. *
  403. ******************************************************************************/
  404.  
  405. static void factor1(n, c, s)
  406. int n, *c, *s;
  407. {
  408.   int i,k;
  409.  
  410.   /* First factor out any 2s. */
  411.  
  412.   n = factor_out(n, 2, c, s);
  413.  
  414.   /* Now any odd factors. */
  415.  
  416.   k = (int)sqrt((DBL)n) + 1;
  417.  
  418.   for (i = 3; (n > 1) && (i <= k); i += 2)
  419.   {
  420.     if (!(n % i))
  421.     {
  422.       n = factor_out(n, i, c, s);
  423.       k = (int)sqrt((DBL)n)+1;
  424.     }
  425.   }
  426.  
  427.   if (n > 1)
  428.   {
  429.     s[(*c)++] = n;
  430.   }
  431. }
  432.  
  433.  
  434.  
  435. /*****************************************************************************
  436. *
  437. * FUNCTION
  438. *
  439. *   binomial
  440. *
  441. * INPUT
  442. *   
  443. * OUTPUT
  444. *   
  445. * RETURNS
  446. *   
  447. * AUTHOR
  448. *
  449. *   Alexander Enzmann
  450. *   
  451. * DESCRIPTION
  452. *
  453. *   Calculate the binomial coefficent of n, r.
  454. *
  455. * CHANGES
  456. *
  457. *   -
  458. *
  459. ******************************************************************************/
  460.  
  461. static long binomial(n, r)
  462. int n, r;
  463. {
  464.   int h,i,j,k,l;
  465.   unsigned long result;
  466.   static int stack1[BINOMSIZE], stack2[BINOMSIZE];
  467.  
  468.   if ((n < 0) || (r < 0) || (r > n))
  469.   {
  470.     result = 0L;
  471.   }
  472.   else
  473.   {
  474.     if (r == n)
  475.     {
  476.       result = 1L;
  477.     }
  478.     else
  479.     {
  480.       if ((r < 15) && (n < 15))
  481.       {
  482.         result = (long)binomials[n][r];
  483.       }
  484.       else
  485.       {
  486.         j = 0;
  487.  
  488.         for (i = r + 1; i <= n; i++)
  489.         {
  490.           stack1[j++] = i;
  491.         }
  492.  
  493.         for (i = 2; i <= (n-r); i++)
  494.         {
  495.           h = 0;
  496.  
  497.           factor1(i, &h, stack2);
  498.  
  499.           for (k = 0; k < h; k++)
  500.           {
  501.             for (l = 0; l < j; l++)
  502.             {
  503.               if (!(stack1[l] % stack2[k]))
  504.               {
  505.                 stack1[l] /= stack2[k];
  506.  
  507.                 goto l1;
  508.               }
  509.             }
  510.           }
  511.  
  512.           /* Error if we get here */
  513. /*        Debug_Info("Failed to factor\n");*/
  514. l1:;
  515.         }
  516.  
  517.         result = 1;
  518.  
  519.         for (i = 0; i < j; i++)
  520.         {
  521.           result *= stack1[i];
  522.         }
  523.       }
  524.     }
  525.   }
  526.  
  527.   return(result);
  528. }
  529.  
  530.  
  531.  
  532. /*****************************************************************************
  533. *
  534. * FUNCTION
  535. *
  536. *   inside
  537. *
  538. * INPUT
  539. *   
  540. * OUTPUT
  541. *   
  542. * RETURNS
  543. *   
  544. * AUTHOR
  545. *
  546. *   Alexander Enzmann
  547. *   
  548. * DESCRIPTION
  549. *
  550. *   -
  551. *
  552. * CHANGES
  553. *
  554. *   -
  555. *
  556. ******************************************************************************/
  557.  
  558. static DBL inside(IPoint, Order, Coeffs)
  559. VECTOR  IPoint;
  560. int Order;
  561. DBL *Coeffs;
  562. {
  563.   DBL x[MAX_ORDER+1], y[MAX_ORDER+1];
  564.   DBL z[MAX_ORDER+1], c, Result;
  565.   int i, j, k, term;
  566.  
  567.   x[0] = 1.0;       y[0] = 1.0;       z[0] = 1.0;
  568.   x[1] = IPoint[X]; y[1] = IPoint[Y]; z[1] = IPoint[Z];
  569.  
  570.   for (i = 2; i <= Order; i++)
  571.   {
  572.     x[i] = x[1] * x[i-1];
  573.     y[i] = y[1] * y[i-1];
  574.     z[i] = z[1] * z[i-1];
  575.   }
  576.  
  577.   Result = 0.0;
  578.  
  579.   term = 0;
  580.  
  581.   for (i = Order; i >= 0; i--)
  582.   {
  583.     for (j=Order-i;j>=0;j--)
  584.     {
  585.       for (k=Order-(i+j);k>=0;k--)
  586.       {
  587.         if ((c = Coeffs[term]) != 0.0)
  588.         {
  589.           Result += c * x[i] * y[j] * z[k];
  590.         }
  591.         term++;
  592.       }
  593.     }
  594.   }
  595.  
  596.   return(Result);
  597. }
  598.  
  599.  
  600.  
  601. /*****************************************************************************
  602. *
  603. * FUNCTION
  604. *
  605. *   intersect
  606. *
  607. * INPUT
  608. *   
  609. * OUTPUT
  610. *   
  611. * RETURNS
  612. *   
  613. * AUTHOR
  614. *
  615. *   Alexander Enzmann
  616. *   
  617. * DESCRIPTION
  618. *
  619. *   Intersection of a ray and an arbitrary polynomial function.
  620. *
  621. * CHANGES
  622. *
  623. *   -
  624. *
  625. ******************************************************************************/
  626.  
  627. static int intersect(ray, Order, Coeffs, Sturm_Flag, Depths)
  628. RAY *ray;
  629. int Order, Sturm_Flag;
  630. DBL *Coeffs, *Depths;
  631. {
  632.   DBL eqn[MAX_ORDER+1];
  633.   DBL t[3][MAX_ORDER+1];
  634.   VECTOR  P, D;
  635.   DBL val;
  636.   int h, i, j, k, i1, j1, k1, term;
  637.   int offset;
  638.  
  639.   /* First we calculate the values of the individual powers
  640.      of x, y, and z as they are represented by the ray */
  641.  
  642.   Assign_Vector(P,ray->Initial);
  643.   Assign_Vector(D,ray->Direction);
  644.  
  645.   for (i = 0; i < 3; i++)
  646.   {
  647.     eqn_v[i][0]  = 1.0;
  648.     eqn_vt[i][0] = 1.0;
  649.   }
  650.  
  651.   eqn_v[0][1] = P[X];
  652.   eqn_v[1][1] = P[Y];
  653.   eqn_v[2][1] = P[Z];
  654.  
  655.   eqn_vt[0][1] = D[X];
  656.   eqn_vt[1][1] = D[Y];
  657.   eqn_vt[2][1] = D[Z];
  658.  
  659.   for (i = 2; i <= Order; i++)
  660.   {
  661.     for (j=0;j<3;j++)
  662.     {
  663.      eqn_v[j][i]  = eqn_v[j][1] * eqn_v[j][i-1];
  664.      eqn_vt[j][i] = eqn_vt[j][1] * eqn_vt[j][i-1];
  665.     }
  666.   }
  667.  
  668.   for (i = 0; i <= Order; i++)
  669.   {
  670.     eqn[i] = 0.0;
  671.   }
  672.  
  673.   /* Now walk through the terms of the polynomial.  As we go
  674.      we substitute the ray equation for each of the variables. */
  675.  
  676.   term = 0;
  677.  
  678.   for (i = Order; i >= 0; i--)
  679.   {
  680.     for (h = 0; h <= i; h++)
  681.     {
  682.       t[0][h] = binomial(i, h) * eqn_vt[0][i-h] * eqn_v[0][h];
  683.     }
  684.  
  685.     for (j = Order-i; j >= 0; j--)
  686.     {
  687.       for (h = 0; h <= j; h++)
  688.       {
  689.         t[1][h] = binomial(j, h) * eqn_vt[1][j-h] * eqn_v[1][h];
  690.       }
  691.  
  692.       for (k = Order-(i+j); k >= 0; k--)
  693.       {
  694.         if (Coeffs[term] != 0)
  695.         {
  696.           for (h = 0; h <= k; h++)
  697.           {
  698.             t[2][h] = binomial(k, h) * eqn_vt[2][k-h] * eqn_v[2][h];
  699.           }
  700.  
  701.           /* Multiply together the three polynomials. */
  702.  
  703.           offset = Order - (i + j + k);
  704.  
  705.           for (i1 = 0; i1 <= i; i1++)
  706.           {
  707.             for (j1=0;j1<=j;j1++)
  708.             {
  709.               for (k1=0;k1<=k;k1++)
  710.               {
  711.                 h = offset + i1 + j1 + k1;
  712.                 val = Coeffs[term];
  713.                 val *= t[0][i1];
  714.                 val *= t[1][j1];
  715.                 val *= t[2][k1];
  716.                 eqn[h] += val;
  717.               }
  718.             }
  719.           }
  720.         }
  721.  
  722.         term++;
  723.       }
  724.     }
  725.   }
  726.  
  727.   for (i = 0, j = Order; i <= Order; i++)
  728.   {
  729.     if (eqn[i] != 0.0)
  730.     {
  731.       break;
  732.     }
  733.     else
  734.     {
  735.       j--;
  736.     }
  737.   }
  738.  
  739.   if (j > 1)
  740.   {
  741.     return(Solve_Polynomial(j, &eqn[i], Depths, Sturm_Flag, ROOT_TOLERANCE));
  742.   }
  743.   else
  744.   {
  745.     return(0);
  746.   }
  747. }
  748.  
  749.  
  750.  
  751. /*****************************************************************************
  752. *
  753. * FUNCTION
  754. *
  755. *   intersect_linear
  756. *
  757. * INPUT
  758. *   
  759. * OUTPUT
  760. *   
  761. * RETURNS
  762. *   
  763. * AUTHOR
  764. *
  765. *   Alexander Enzmann
  766. *   
  767. * DESCRIPTION
  768. *
  769. *   Intersection of a ray and a quadratic.
  770. *
  771. * CHANGES
  772. *
  773. *   -
  774. *
  775. ******************************************************************************/
  776.  
  777. static int intersect_linear(ray, Coeffs, Depths)
  778. RAY *ray;
  779. DBL *Coeffs, *Depths;
  780. {
  781.   DBL t0, t1, *a = Coeffs;
  782.  
  783.   t0 = a[0] * ray->Initial[X] + a[1] * ray->Initial[Y] + a[2] * ray->Initial[Z];
  784.   t1 = a[0] * ray->Direction[X] + a[1] * ray->Direction[Y] +
  785.  
  786.   a[2] * ray->Direction[Z];
  787.  
  788.   if (fabs(t1) < EPSILON)
  789.   {
  790.     return(0);
  791.   }
  792.  
  793.   Depths[0] = -(a[3] + t0) / t1;
  794.  
  795.   return(1);
  796. }
  797.  
  798.  
  799.  
  800. /*****************************************************************************
  801. *
  802. * FUNCTION
  803. *
  804. *   intersect_quadratic
  805. *
  806. * INPUT
  807. *   
  808. * OUTPUT
  809. *   
  810. * RETURNS
  811. *   
  812. * AUTHOR
  813. *
  814. *   Alexander Enzmann
  815. *   
  816. * DESCRIPTION
  817. *
  818. *   Intersection of a ray and a quadratic.
  819. *
  820. * CHANGES
  821. *
  822. *   -
  823. *
  824. ******************************************************************************/
  825.  
  826. static int intersect_quadratic(ray, Coeffs, Depths)
  827. RAY *ray;
  828. DBL *Coeffs, *Depths;
  829. {
  830.   DBL x, y, z, x2, y2, z2;
  831.   DBL xx, yy, zz, xx2, yy2, zz2;
  832.   DBL *a, ac, bc, cc, d, t;
  833.  
  834.   x  = ray->Initial[X];
  835.   y  = ray->Initial[Y];
  836.   z  = ray->Initial[Z];
  837.  
  838.   xx = ray->Direction[X];
  839.   yy = ray->Direction[Y];
  840.   zz = ray->Direction[Z];
  841.  
  842.   x2  = x * x;    y2  = y * y;    z2 = z * z;
  843.   xx2 = xx * xx;  yy2 = yy * yy;  zz2 = zz * zz;
  844.  
  845.   a = Coeffs;
  846.  
  847.   /*
  848.    * Determine the coefficients of t^n, where the line is represented
  849.    * as (x,y,z) + (xx,yy,zz)*t.
  850.    */
  851.  
  852.   ac = (a[0]*xx2 + a[1]*xx*yy + a[2]*xx*zz + a[4]*yy2 + a[5]*yy*zz + a[7]*zz2);
  853.  
  854.   bc = (2*a[0]*x*xx + a[1]*(x*yy + xx*y) + a[2]*(x*zz + xx*z) +
  855.         a[3]*xx + 2*a[4]*y*yy + a[5]*(y*zz + yy*z) + a[6]*yy +
  856.         2*a[7]*z*zz + a[8]*zz);
  857.  
  858.   cc = a[0]*x2 + a[1]*x*y + a[2]*x*z + a[3]*x + a[4]*y2 +
  859.        a[5]*y*z + a[6]*y + a[7]*z2 + a[8]*z + a[9];
  860.  
  861.   if (fabs(ac) < COEFF_LIMIT)
  862.   {
  863.     if (fabs(bc) < COEFF_LIMIT)
  864.     {
  865.       return(0);
  866.     }
  867.  
  868.     Depths[0] = -cc / bc;
  869.  
  870.     return(1);
  871.   }
  872.  
  873.   /*
  874.    * Solve the quadratic formula & return results that are
  875.    * within the correct interval.
  876.    */
  877.  
  878.   d = bc * bc - 4.0 * ac * cc;
  879.  
  880.   if (d < 0.0)
  881.   {
  882.     return(0);
  883.   }
  884.  
  885.   d = sqrt(d);
  886.  
  887.   bc = -bc;
  888.  
  889.   t = (1.0 / (2.0 * ac));
  890.   Depths[0] = (bc + d) * t;
  891.   Depths[1] = (bc - d) * t;
  892.  
  893. /*
  894.   t = 2.0 * ac;
  895.   Depths[0] = (bc + d) / t;
  896.   Depths[1] = (bc - d) / t;
  897. */
  898.  
  899.   return(2);
  900. }
  901.  
  902.  
  903.  
  904. /*****************************************************************************
  905. *
  906. * FUNCTION
  907. *
  908. *   normal0
  909. *
  910. * INPUT
  911. *   
  912. * OUTPUT
  913. *   
  914. * RETURNS
  915. *   
  916. * AUTHOR
  917. *
  918. *   Alexander Enzmann
  919. *   
  920. * DESCRIPTION
  921. *
  922. *   Normal to a polynomial (used for polynomials with order > 4).
  923. *
  924. * CHANGES
  925. *
  926. *   -
  927. *
  928. ******************************************************************************/
  929.  
  930. static void normal0(Result, Order, Coeffs, IPoint)
  931. VECTOR Result;
  932. int Order;
  933. DBL *Coeffs;
  934. VECTOR IPoint;
  935.  {
  936.   int i, j, k, term;
  937.   DBL val, *a, x[MAX_ORDER+1], y[MAX_ORDER+1], z[MAX_ORDER+1];
  938.  
  939.   x[0] = 1.0;
  940.   y[0] = 1.0;
  941.   z[0] = 1.0;
  942.  
  943.   x[1] = IPoint[X];
  944.   y[1] = IPoint[Y];
  945.   z[1] = IPoint[Z];
  946.  
  947.   for (i = 2; i <= Order; i++)
  948.   {
  949.     x[i] = IPoint[X] * x[i-1];
  950.     y[i] = IPoint[Y] * y[i-1];
  951.     z[i] = IPoint[Z] * z[i-1];
  952.   }
  953.  
  954.   a = Coeffs;
  955.  
  956.   term = 0;
  957.  
  958.   Make_Vector(Result, 0.0, 0.0, 0.0);
  959.  
  960.   for (i = Order; i >= 0; i--)
  961.   {
  962.     for (j = Order-i; j >= 0; j--)
  963.     {
  964.       for (k = Order-(i+j); k >= 0; k--)
  965.       {
  966.         if (i >= 1)
  967.         {
  968.           val = x[i-1] * y[j] * z[k];
  969.           Result[X] += i * a[term] * val;
  970.         }
  971.  
  972.         if (j >= 1)
  973.         {
  974.           val = x[i] * y[j-1] * z[k];
  975.           Result[Y] += j * a[term] * val;
  976.         }
  977.  
  978.         if (k >= 1)
  979.         {
  980.           val = x[i] * y[j] * z[k-1];
  981.           Result[Z] += k * a[term] * val;
  982.         }
  983.  
  984.         term++;
  985.       }
  986.     }
  987.   }
  988. }
  989.  
  990.  
  991.  
  992. /*****************************************************************************
  993. *
  994. * FUNCTION
  995. *
  996. *   nromal1
  997. *
  998. * INPUT
  999. *   
  1000. * OUTPUT
  1001. *   
  1002. * RETURNS
  1003. *   
  1004. * AUTHOR
  1005. *
  1006. *   Alexander Enzmann
  1007. *   
  1008. * DESCRIPTION
  1009. *
  1010. *   Normal to a polynomial (for polynomials of order <= 4).
  1011. *
  1012. * CHANGES
  1013. *
  1014. *   -
  1015. *
  1016. ******************************************************************************/
  1017.  
  1018. static void normal1(Result, Order, Coeffs, IPoint)
  1019. VECTOR Result;
  1020. int Order;
  1021. DBL *Coeffs;
  1022. VECTOR IPoint;
  1023. {
  1024.   DBL *a, x, y, z, x2, y2, z2, x3, y3, z3;
  1025.  
  1026.   a = Coeffs;
  1027.  
  1028.   x = IPoint[X];
  1029.   y = IPoint[Y];
  1030.   z = IPoint[Z];
  1031.  
  1032.   switch (Order)
  1033.   {
  1034.     case 1:
  1035.  
  1036.       /* Linear partial derivatives */
  1037.  
  1038.       Make_Vector(Result, a[0], a[1], a[2])
  1039.  
  1040.       break;
  1041.  
  1042.     case 2:
  1043.  
  1044.       /* Quadratic partial derivatives */
  1045.  
  1046.       Result[X] = 2*a[0]*x+a[1]*y+a[2]*z+a[3];
  1047.       Result[Y] = a[1]*x+2*a[4]*y+a[5]*z+a[6];
  1048.       Result[Z] = a[2]*x+a[5]*y+2*a[7]*z+a[8];
  1049.  
  1050.       break;
  1051.  
  1052.     case 3:
  1053.  
  1054.         x2 = x * x;  y2 = y * y;  z2 = z * z;
  1055.  
  1056.         /* Cubic partial derivatives */
  1057.  
  1058.         Result[X] = 3*a[0]*x2 + 2*x*(a[1]*y + a[2]*z + a[3]) + a[4]*y2 +
  1059.                     y*(a[5]*z + a[6]) + a[7]*z2 + a[8]*z + a[9];
  1060.         Result[Y] = a[1]*x2 + x*(2*a[4]*y + a[5]*z + a[6]) + 3*a[10]*y2 +
  1061.                     2*y*(a[11]*z + a[12]) + a[13]*z2 + a[14]*z + a[15];
  1062.         Result[Z] = a[2]*x2 + x*(a[5]*y + 2*a[7]*z + a[8]) + a[11]*y2 +
  1063.                     y*(2*a[13]*z + a[14]) + 3*a[16]*z2 + 2*a[17]*z + a[18];
  1064.  
  1065.         break;
  1066.  
  1067.     case 4:
  1068.  
  1069.       /* Quartic partial derivatives */
  1070.  
  1071.       x2 = x * x;  y2 = y * y;  z2 = z * z;
  1072.       x3 = x * x2; y3 = y * y2; z3 = z * z2;
  1073.  
  1074.       Result[X] = 4*a[ 0]*x3+3*x2*(a[ 1]*y+a[ 2]*z+a[ 3])+
  1075.                   2*x*(a[ 4]*y2+y*(a[ 5]*z+a[ 6])+a[ 7]*z2+a[ 8]*z+a[ 9])+
  1076.                   a[10]*y3+y2*(a[11]*z+a[12])+y*(a[13]*z2+a[14]*z+a[15])+
  1077.                   a[16]*z3+a[17]*z2+a[18]*z+a[19];
  1078.       Result[Y] = a[ 1]*x3+x2*(2*a[ 4]*y+a[ 5]*z+a[ 6])+
  1079.                   x*(3*a[10]*y2+2*y*(a[11]*z+a[12])+a[13]*z2+a[14]*z+a[15])+
  1080.                   4*a[20]*y3+3*y2*(a[21]*z+a[22])+2*y*(a[23]*z2+a[24]*z+a[25])+
  1081.                   a[26]*z3+a[27]*z2+a[28]*z+a[29];
  1082.       Result[Z] = a[ 2]*x3+x2*(a[ 5]*y+2*a[ 7]*z+a[ 8])+
  1083.                   x*(a[11]*y2+y*(2*a[13]*z+a[14])+3*a[16]*z2+2*a[17]*z+a[18])+
  1084.                   a[21]*y3+y2*(2*a[23]*z+a[24])+y*(3*a[26]*z2+2*a[27]*z+a[28])+
  1085.                   4*a[30]*z3+3*a[31]*z2+2*a[32]*z+a[33];
  1086.   }
  1087. }
  1088.  
  1089.  
  1090.  
  1091. /*****************************************************************************
  1092. *
  1093. * FUNCTION
  1094. *
  1095. *   Inside_Poly
  1096. *
  1097. * INPUT
  1098. *   
  1099. * OUTPUT
  1100. *   
  1101. * RETURNS
  1102. *   
  1103. * AUTHOR
  1104. *
  1105. *   Alexander Enzmann
  1106. *
  1107. * DESCRIPTION
  1108. *
  1109. *   -
  1110. *
  1111. * CHANGES
  1112. *
  1113. *   -
  1114. *
  1115. ******************************************************************************/
  1116.  
  1117. static int Inside_Poly (IPoint, Object)
  1118. VECTOR IPoint;
  1119. OBJECT *Object;
  1120. {
  1121.   VECTOR  New_Point;
  1122.   DBL Result;
  1123.  
  1124.   /* Transform the point into polynomial's space */
  1125.  
  1126.   MInvTransPoint(New_Point, IPoint, ((POLY *)Object)->Trans);
  1127.  
  1128.   Result = inside(New_Point, ((POLY *)Object)->Order, ((POLY *)Object)->Coeffs);
  1129.  
  1130.   if (Result < INSIDE_TOLERANCE)
  1131.   {
  1132.     return ((int)(!Test_Flag(Object, INVERTED_FLAG)));
  1133.   }
  1134.   else
  1135.   {
  1136.     return ((int)Test_Flag(Object, INVERTED_FLAG));
  1137.   }
  1138. }
  1139.  
  1140.  
  1141.  
  1142. /*****************************************************************************
  1143. *
  1144. * FUNCTION
  1145. *
  1146. *   Poly_Normal
  1147. *
  1148. * INPUT
  1149. *   
  1150. * OUTPUT
  1151. *
  1152. * RETURNS
  1153. *   
  1154. * AUTHOR
  1155. *
  1156. *   Alexander Enzmann
  1157. *   
  1158. * DESCRIPTION
  1159. *
  1160. *   Normal to a polynomial.
  1161. *
  1162. * CHANGES
  1163. *
  1164. *   -
  1165. *
  1166. ******************************************************************************/
  1167.  
  1168. static void Poly_Normal(Result, Object, Inter)
  1169. OBJECT *Object;
  1170. VECTOR Result;
  1171. INTERSECTION *Inter;
  1172. {
  1173.   DBL val;
  1174.   VECTOR  New_Point;
  1175.   POLY *Shape = (POLY *)Object;
  1176.  
  1177.   /* Transform the point into the polynomials space. */
  1178.  
  1179.   MInvTransPoint(New_Point, Inter->IPoint, Shape->Trans);
  1180.  
  1181.   if (((POLY *)Object)->Order > 4)
  1182.   {
  1183.     normal0(Result, Shape->Order, Shape->Coeffs, New_Point);
  1184.   }
  1185.   else
  1186.   {
  1187.     normal1(Result, Shape->Order, Shape->Coeffs, New_Point);
  1188.   }
  1189.  
  1190.   /* Transform back to world space. */
  1191.  
  1192.   MTransNormal(Result, Result, Shape->Trans);
  1193.  
  1194.   /* Normalize (accounting for the possibility of a 0 length normal). */
  1195.  
  1196.   VDot(val, Result, Result);
  1197.  
  1198.   if (val > 0.0)
  1199.   {
  1200.     val = 1.0 / sqrt(val);
  1201.  
  1202.     VScaleEq(Result, val);
  1203.   }
  1204.   else
  1205.   {
  1206.     Make_Vector(Result, 1.0, 0.0, 0.0)
  1207.   }
  1208. }
  1209.  
  1210.  
  1211.  
  1212. /*****************************************************************************
  1213. *
  1214. * FUNCTION
  1215. *
  1216. *   Translate_Poly
  1217. *
  1218. * INPUT
  1219. *   
  1220. * OUTPUT
  1221. *   
  1222. * RETURNS
  1223. *   
  1224. * AUTHOR
  1225. *
  1226. *   Alexander Enzmann
  1227. *   
  1228. * DESCRIPTION
  1229. *
  1230. *   -
  1231. *
  1232. * CHANGES
  1233. *
  1234. *   -
  1235. *
  1236. ******************************************************************************/
  1237.  
  1238. static void Translate_Poly (Object, Vector, Trans)
  1239. OBJECT *Object;
  1240. VECTOR Vector;
  1241. TRANSFORM *Trans;
  1242. {
  1243.   Transform_Poly(Object, Trans);
  1244. }
  1245.  
  1246.  
  1247.  
  1248. /*****************************************************************************
  1249. *
  1250. * FUNCTION
  1251. *
  1252. *   Rotate_Poly
  1253. *
  1254. * INPUT
  1255. *   
  1256. * OUTPUT
  1257. *   
  1258. * RETURNS
  1259. *   
  1260. * AUTHOR
  1261. *
  1262. *   Alexander Enzmann
  1263. *   
  1264. * DESCRIPTION
  1265. *
  1266. *   -
  1267. *
  1268. * CHANGES
  1269. *
  1270. *   -
  1271. *
  1272. ******************************************************************************/
  1273.  
  1274. static void Rotate_Poly (Object, Vector, Trans)
  1275. OBJECT *Object;
  1276. VECTOR Vector;
  1277. TRANSFORM *Trans;
  1278. {
  1279.   Transform_Poly(Object, Trans);
  1280. }
  1281.  
  1282.  
  1283.  
  1284. /*****************************************************************************
  1285. *
  1286. * FUNCTION
  1287. *
  1288. *   Scale_Poly
  1289. *
  1290. * INPUT
  1291. *   
  1292. * OUTPUT
  1293. *   
  1294. * RETURNS
  1295. *   
  1296. * AUTHOR
  1297. *
  1298. *   Alexander Enzmann
  1299. *   
  1300. * DESCRIPTION
  1301. *
  1302. *   -
  1303. *
  1304. * CHANGES
  1305. *
  1306. *   -
  1307. *
  1308. ******************************************************************************/
  1309.  
  1310. static void Scale_Poly (Object, Vector, Trans)
  1311. OBJECT *Object;
  1312. VECTOR Vector;
  1313. TRANSFORM *Trans;
  1314. {
  1315.   Transform_Poly(Object, Trans);
  1316. }
  1317.  
  1318.  
  1319.  
  1320. /*****************************************************************************
  1321. *
  1322. * FUNCTION
  1323. *
  1324. *   Transform_Poly
  1325. *
  1326. * INPUT
  1327. *   
  1328. * OUTPUT
  1329. *   
  1330. * RETURNS
  1331. *   
  1332. * AUTHOR
  1333. *
  1334. *   Alexander Enzmann
  1335. *   
  1336. * DESCRIPTION
  1337. *
  1338. *   -
  1339. *
  1340. * CHANGES
  1341. *
  1342. *   -
  1343. *
  1344. ******************************************************************************/
  1345.  
  1346. static void Transform_Poly(Object,Trans)
  1347. OBJECT *Object;
  1348. TRANSFORM *Trans;
  1349. {
  1350.   Compose_Transforms(((POLY *)Object)->Trans, Trans);
  1351.  
  1352.   Compute_Poly_BBox((POLY *)Object);
  1353. }
  1354.  
  1355.  
  1356.  
  1357. /*****************************************************************************
  1358. *
  1359. * FUNCTION
  1360. *
  1361. *   Invert_Poly
  1362. *
  1363. * INPUT
  1364. *   
  1365. * OUTPUT
  1366. *   
  1367. * RETURNS
  1368. *   
  1369. * AUTHOR
  1370. *
  1371. *   Alexander Enzmann
  1372. *   
  1373. * DESCRIPTION
  1374. *
  1375. *   -
  1376. *
  1377. * CHANGES
  1378. *
  1379. *   -
  1380. *
  1381. ******************************************************************************/
  1382.  
  1383. static void Invert_Poly(Object)
  1384. OBJECT *Object;
  1385. {
  1386.   Invert_Flag(Object, INVERTED_FLAG);
  1387. }
  1388.  
  1389.  
  1390.  
  1391. /*****************************************************************************
  1392. *
  1393. * FUNCTION
  1394. *
  1395. *   Create_Poly
  1396. *
  1397. * INPUT
  1398. *   
  1399. * OUTPUT
  1400. *   
  1401. * RETURNS
  1402. *   
  1403. * AUTHOR
  1404. *
  1405. *   Alexander Enzmann
  1406. *   
  1407. * DESCRIPTION
  1408. *
  1409. *   -
  1410. *
  1411. * CHANGES
  1412. *
  1413. *   -
  1414. *
  1415. ******************************************************************************/
  1416.  
  1417. POLY *Create_Poly(Order)
  1418. int Order;
  1419. {
  1420.   POLY *New;
  1421.   int i;
  1422.  
  1423.   New = (POLY *)POV_MALLOC(sizeof (POLY), "poly");
  1424.  
  1425.   INIT_OBJECT_FIELDS(New,POLY_OBJECT, &Poly_Methods);
  1426.  
  1427.   New->Order = Order;
  1428.  
  1429.   New->Trans = Create_Transform();
  1430.  
  1431.   New->Coeffs = (DBL *)POV_MALLOC(term_counts(Order) * sizeof(DBL), "coefficients for POLY");
  1432.  
  1433.   for (i = 0; i < term_counts(Order); i++)
  1434.   {
  1435.     New->Coeffs[i] = 0.0;
  1436.   }
  1437.  
  1438.   return(New);
  1439. }
  1440.  
  1441.  
  1442.  
  1443. /*****************************************************************************
  1444. *
  1445. * FUNCTION
  1446. *
  1447. *   Copy_Poly
  1448. *
  1449. * INPUT
  1450. *   
  1451. * OUTPUT
  1452. *   
  1453. * RETURNS
  1454. *   
  1455. * AUTHOR
  1456. *
  1457. *   Alexander Enzmann
  1458. *   
  1459. * DESCRIPTION
  1460. *
  1461. *   Make a copy of a polynomial object.
  1462. *
  1463. * CHANGES
  1464. *
  1465. *   -
  1466. *
  1467. ******************************************************************************/
  1468.  
  1469. static void *Copy_Poly(Object)
  1470. OBJECT *Object;
  1471. {
  1472.   POLY *New;
  1473.   int i;
  1474.  
  1475.   New = Create_Poly (((POLY *)Object)->Order);
  1476.  
  1477.   /* Get rid of transform created in Create_Poly. */
  1478.  
  1479.   Destroy_Transform(New->Trans);
  1480.  
  1481.   Copy_Flag(New, Object, STURM_FLAG);
  1482.   Copy_Flag(New, Object, INVERTED_FLAG);
  1483.  
  1484.   New->Trans = Copy_Transform(((POLY *)Object)->Trans);
  1485.  
  1486.   for (i = 0; i < term_counts(New->Order); i++)
  1487.   {
  1488.     New->Coeffs[i] = ((POLY *)Object)->Coeffs[i];
  1489.   }
  1490.  
  1491.   return (New);
  1492. }
  1493.  
  1494.  
  1495.  
  1496. /*****************************************************************************
  1497. *
  1498. * FUNCTION
  1499. *
  1500. *   Destroy_Poly
  1501. *
  1502. * INPUT
  1503. *   
  1504. * OUTPUT
  1505. *   
  1506. * RETURNS
  1507. *   
  1508. * AUTHOR
  1509. *
  1510. *   Alexander Enzmann
  1511. *   
  1512. * DESCRIPTION
  1513. *
  1514. *   -
  1515. *
  1516. * CHANGES
  1517. *
  1518. *   -
  1519. *
  1520. ******************************************************************************/
  1521.  
  1522. static void Destroy_Poly(Object)
  1523. OBJECT *Object;
  1524. {
  1525.   Destroy_Transform (((POLY *)Object)->Trans);
  1526.  
  1527.   POV_FREE (((POLY *)Object)->Coeffs);
  1528.  
  1529.   POV_FREE (Object);
  1530. }
  1531.  
  1532.  
  1533.  
  1534. /*****************************************************************************
  1535. *
  1536. * FUNCTION
  1537. *
  1538. *   Compute_Poly_BBox
  1539. *
  1540. * INPUT
  1541. *
  1542. *   Poly - Poly
  1543. *   
  1544. * OUTPUT
  1545. *
  1546. *   Poly
  1547. *
  1548. * RETURNS
  1549. *   
  1550. * AUTHOR
  1551. *
  1552. *   Dieter Bayer
  1553. *   
  1554. * DESCRIPTION
  1555. *
  1556. *   Calculate the bounding box of a poly.
  1557. *
  1558. * CHANGES
  1559. *
  1560. *   Aug 1994 : Creation.
  1561. *
  1562. ******************************************************************************/
  1563.  
  1564. void Compute_Poly_BBox(Poly)
  1565. POLY *Poly;
  1566. {
  1567.   Make_BBox(Poly->BBox, -BOUND_HUGE/2, -BOUND_HUGE/2, -BOUND_HUGE/2, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
  1568.  
  1569.   if (Poly->Clip != NULL)
  1570.   {
  1571.     Poly->BBox = Poly->Clip->BBox;
  1572.   }
  1573. }
  1574.  
  1575.